Reciprocal Inverse Gaussian Distribution (recipinvgauss)

Reciprocal Inverse Gaussian Distribution (recipinvgauss)#

The reciprocal inverse Gaussian distribution is the distribution of a positive random variable whose reciprocal is inverse Gaussian. It naturally appears when you model times with an inverse Gaussian and then switch to modeling the corresponding rates (reciprocals of times).

What you’ll learn#

  • what recipinvgauss models and when it’s appropriate

  • the PDF/CDF (with LaTeX) and its relationship to invgauss

  • key moments (mean/variance/skew/kurtosis), MGF/CF, and entropy

  • how the parameter mu changes the shape

  • core derivations: PDF via change-of-variables, moments via the MGF, likelihood + MLE

  • NumPy-only sampling (via an inverse-Gaussian sampler) + Monte Carlo validation

  • SciPy usage: scipy.stats.recipinvgauss (pdf, cdf, rvs, fit)

  • hypothesis testing, Bayesian modeling patterns, and generative modeling ideas

import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio

from scipy import special, stats

# Plotly rendering (CKC convention)
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

# Reproducibility
rng = np.random.default_rng(7)

np.set_printoptions(precision=4, suppress=True)

1) Title & Classification#

  • Name: recipinvgauss (reciprocal inverse Gaussian)

  • Type: Continuous

  • Support (standard form): (x > 0)

  • Parameter space (SciPy):

    • shape (\mu > 0)

    • loc (\in \mathbb{R})

    • scale (> 0)

We write (standardized form):

[ X \sim \mathrm{RIG}(\mu)\qquad (\mu>0), ]

and SciPy’s location–scale version:

[ X = \mathrm{loc} + \mathrm{scale},Y,\quad Y \sim \mathrm{RIG}(\mu). ]

2) Intuition & Motivation#

2.1 What it models#

A clean way to understand recipinvgauss is through a transformation:

  • Let (Y) be inverse Gaussian.

  • Define (X = 1/Y).

Then (X) has the reciprocal inverse Gaussian distribution.

This matters because many problems naturally flip between a quantity and its reciprocal:

  • time (\leftrightarrow) rate

  • variance (\leftrightarrow) precision

  • speed (\leftrightarrow) travel time per unit distance

If inverse Gaussian is a good model for a positive time variable, recipinvgauss becomes a principled model for the corresponding rate.

2.2 Real-world use cases (examples)#

  • Reliability / survival / first-passage times: inverse Gaussian is a classic first-passage-time model; reciprocals correspond to modeling rates or intensities derived from those times.

  • Hierarchical Bayesian models: if you model a positive scale parameter with an inverse Gaussian, then the corresponding precision parameter follows a reciprocal inverse Gaussian.

  • Variance–mean mixtures: recipinvgauss is a special case of the generalized inverse Gaussian family, which is widely used as a mixing distribution to build heavy-tailed marginals.

2.3 Relations to other distributions#

  • Inverse Gaussian: if (X \sim \mathrm{RIG}(\mu)), then (1/X) is inverse Gaussian.

  • Generalized inverse Gaussian (GIG): recipinvgauss is a (\mathrm{GIG}(p=1/2,,a=1,,b=1/\mu^2)) distribution.

3) Formal Definition#

SciPy’s standardized recipinvgauss(mu) has PDF:

[ f(x\mid\mu) = \frac{1}{\sqrt{2\pi x}}\exp\left(-\frac{(1-\mu x)^2}{2\mu^2 x}\right),\qquad x>0,;\mu>0. ]

A convenient CDF expression uses the standard normal CDF (\Phi):

[ F(x\mid\mu) = \Phi!\left(\frac{\mu x - 1}{\mu\sqrt{x}}\right) ; -; \exp!\left(\frac{2}{\mu}\right) ,\Phi!\left(-\frac{\mu x + 1}{\mu\sqrt{x}}\right),\qquad x>0. ]

Relationship to invgauss#

If (Y \sim \mathrm{InvGauss}(\mu)) in SciPy’s standardized form and (X = 1/Y), then:

[ X \sim \mathrm{RIG}(\mu), \qquad F_X(x) = \mathbb{P}(X\le x) = \mathbb{P}(Y\ge 1/x) = 1 - F_Y(1/x). ]

def recipinvgauss_logpdf(x: np.ndarray, mu: float) -> np.ndarray:
    '''Log-PDF of the standardized reciprocal inverse Gaussian.

    Parameters
    ----------
    x:
        points where to evaluate (array-like)
    mu:
        shape parameter, mu > 0
    '''

    x = np.asarray(x, dtype=float)
    mu = float(mu)
    if mu <= 0:
        raise ValueError('mu must be > 0')

    out = np.full_like(x, -np.inf, dtype=float)
    mask = x > 0

    xm = x[mask]
    out[mask] = (
        -0.5 * np.log(2.0 * np.pi)
        - 0.5 * np.log(xm)
        - ((1.0 - mu * xm) ** 2) / (2.0 * mu**2 * xm)
    )
    return out


def recipinvgauss_pdf(x: np.ndarray, mu: float) -> np.ndarray:
    return np.exp(recipinvgauss_logpdf(x, mu))


def recipinvgauss_cdf(x: np.ndarray, mu: float) -> np.ndarray:
    '''CDF via the closed-form normal-CDF expression.

    Uses log-space for the exp(2/mu)*Phi(.) term for stability.
    '''

    x = np.asarray(x, dtype=float)
    mu = float(mu)
    if mu <= 0:
        raise ValueError('mu must be > 0')

    out = np.zeros_like(x, dtype=float)
    mask = x > 0
    xm = x[mask]

    a = (mu * xm - 1.0) / (mu * np.sqrt(xm))
    b = -(mu * xm + 1.0) / (mu * np.sqrt(xm))

    term1 = special.ndtr(a)
    term2 = np.exp(2.0 / mu + special.log_ndtr(b))

    out[mask] = term1 - term2
    return np.clip(out, 0.0, 1.0)

# Quick sanity checks vs SciPy
mu0 = 2.5
x_grid = np.logspace(-3, 2, 200)

pdf_max_err = np.max(
    np.abs(recipinvgauss_pdf(x_grid, mu0) - stats.recipinvgauss(mu0).pdf(x_grid))
)

cdf_max_err = np.max(
    np.abs(recipinvgauss_cdf(x_grid, mu0) - stats.recipinvgauss(mu0).cdf(x_grid))
)

pdf_max_err, cdf_max_err

(1.1102230246251565e-16, 3.3306690738754696e-16)

4) Moments & Properties#

4.1 Mean, variance, skewness, kurtosis#

For (X \sim \mathrm{RIG}(\mu)):

  • Mean [ \mathbb{E}[X] = 1 + \frac{1}{\mu}. ]

  • Variance [ \mathrm{Var}(X) = 2 + \frac{1}{\mu}. ]

  • Skewness [ \gamma_1 = \frac{8 + 3/\mu}{\left(2 + 1/\mu\right)^{3/2}}. ]

  • Excess kurtosis (kurtosis minus 3) [ \gamma_2 = \frac{3\left(16 + 5/\mu\right)}{\left(2 + 1/\mu\right)^2}. ]

A useful extra property is the mode (maximizer of the PDF):

[ \mathrm{mode}(X) = \frac{\sqrt{1+4/\mu^2}-1}{2}. ]

4.2 MGF and characteristic function#

The MGF exists for (t < 1/2) and has a simple closed form:

[ M_X(t) = \mathbb{E}[e^{tX}] = (1-2t)^{-1/2},\exp\left(\frac{1-\sqrt{1-2t}}{\mu}\right),\qquad t < 1/2. ]

The characteristic function follows by substituting (t \mapsto i t):

[ \varphi_X(t) = (1-2 i t)^{-1/2},\exp\left(\frac{1-\sqrt{1-2 i t}}{\mu}\right). ]

4.3 Raw moments via Bessel functions (GIG view)#

Because recipinvgauss is a special case of the GIG family, its raw moments can be written with the modified Bessel (K) function:

[ \mathbb{E}[X^r] = \mu^{-r},\frac{K_{r+1/2}(1/\mu)}{K_{1/2}(1/\mu)}. ]

For half-integer orders, (K_{n+1/2}) reduces to (e^{-z}) times a polynomial in (1/z), which is why the low-order moments above become simple rational expressions in (\mu).

4.4 Entropy#

The differential entropy is

[ H(X) = -\int_0^\infty f(x),\log f(x),dx. ]

For this distribution SciPy computes entropy by numerical integration (.entropy()), and that’s typically the most practical approach.

def recipinvgauss_moments(mu: float) -> dict:
    mu = float(mu)
    if mu <= 0:
        raise ValueError('mu must be > 0')

    mean = 1.0 + 1.0 / mu
    var = 2.0 + 1.0 / mu
    skew = (8.0 + 3.0 / mu) / (var ** 1.5)
    excess_kurt = 3.0 * (16.0 + 5.0 / mu) / (var**2)
    mode = 0.5 * (np.sqrt(1.0 + 4.0 / mu**2) - 1.0)

    def mgf(t: float) -> float:
        t = float(t)
        if t >= 0.5:
            raise ValueError('MGF exists only for t < 1/2')
        return (1.0 - 2.0 * t) ** (-0.5) * np.exp((1.0 - np.sqrt(1.0 - 2.0 * t)) / mu)

    def cf(t: float) -> complex:
        t = float(t)
        z = 1.0 - 2.0j * t
        return z ** (-0.5) * np.exp((1.0 - np.sqrt(z)) / mu)

    return {
        'mean': mean,
        'var': var,
        'skew': skew,
        'excess_kurt': excess_kurt,
        'mode': mode,
        'mgf': mgf,
        'cf': cf,
    }


m = recipinvgauss_moments(mu0)
scipy_stats = stats.recipinvgauss.stats(mu0, moments='mvsk')

m['mean'], scipy_stats[0], m['var'], scipy_stats[1], m['skew'], scipy_stats[2], m['excess_kurt'], scipy_stats[3]

(1.4,
 1.4000000000102597,
 2.4,
 2.399999998838463,
 2.4744060267436274,
 2.4744060304700746,
 9.375,
 9.374999944645813)

# Monte Carlo check (mean/var + MGF at a few t)
mu0 = 2.5
n = 300_000
samples = stats.recipinvgauss(mu0).rvs(size=n, random_state=rng)

mc_mean = samples.mean()
mc_var = samples.var(ddof=0)

for t in [0.1, 0.2, -0.2]:
    mc = np.mean(np.exp(t * samples))
    th = recipinvgauss_moments(mu0)['mgf'](t)
    print(f't={t:+.2f}  MC={mc:.6f}  theory={th:.6f}')

mc_mean, recipinvgauss_moments(mu0)['mean'], mc_var, recipinvgauss_moments(mu0)['var']

t=+0.10  MC=1.166121  theory=1.166259
t=+0.20  MC=1.412576  theory=1.412801
t=-0.20  MC=0.785599  theory=0.785431

(1.398831941273997, 1.4, 2.399184268895933, 2.4)

# Entropy (SciPy computes it numerically)
for mu in [0.5, 1.0, 2.5]:
    h = stats.recipinvgauss(mu).entropy()
    print(f'mu={mu:g}  entropy={h:.6f}')

mu=0.5  entropy=1.868685
mu=1  entropy=1.599603
mu=2.5  entropy=1.306416

5) Parameter Interpretation#

A handy interpretation comes from the reciprocal relationship:

[ X \sim \mathrm{RIG}(\mu) \quad\Longleftrightarrow\quad \frac{1}{X} \sim \mathrm{InvGauss}(\mu). ]

So (\mu) is the mean of the reciprocal:

[ \mathbb{E}\left[\frac{1}{X}\right] = \mu. ]

5.1 Shape changes#

  • Smaller (\mu) pushes mass to larger values (heavier right tail), because the reciprocal has smaller mean.

  • Larger (\mu) shifts the distribution toward smaller typical values (mode shrinks like (\approx 1/\mu^2)), while the mean approaches 1: (\mathbb{E}[X] \to 1) and (\mathrm{Var}(X) \to 2) as (\mu\to\infty).

5.2 Location–scale parameters (SciPy)#

SciPy’s loc and scale transform the variable:

[ X = \mathrm{loc} + \mathrm{scale},Y,\quad Y \sim \mathrm{RIG}(\mu). ]

  • scale stretches/shrinks the distribution (mean and std scale accordingly)

  • loc shifts the support to (x > \mathrm{loc})

# PDF shape for different mu
mus = [0.3, 0.7, 1.5, 4.0]

x = np.linspace(1e-4, 12, 700)
fig = go.Figure()

for mu in mus:
    fig.add_trace(go.Scatter(x=x, y=recipinvgauss_pdf(x, mu), mode='lines', name=f'mu={mu:g}'))

fig.update_layout(
    title='recipinvgauss PDF for various μ',
    xaxis_title='x',
    yaxis_title='f(x)',
)
fig.show()

# Mean and mode as functions of mu
mu_grid = np.linspace(0.2, 6.0, 200)
means = np.array([recipinvgauss_moments(mu)['mean'] for mu in mu_grid])
modes = np.array([recipinvgauss_moments(mu)['mode'] for mu in mu_grid])

fig = go.Figure()
fig.add_trace(go.Scatter(x=mu_grid, y=means, mode='lines', name='mean'))
fig.add_trace(go.Scatter(x=mu_grid, y=modes, mode='lines', name='mode'))
fig.update_layout(
    title='Mean and mode vs μ',
    xaxis_title='μ',
    yaxis_title='value',
)
fig.show()

6) Derivations#

6.1 Expectation#

From the MGF,

[ M(t) = (1-2t)^{-1/2}\exp\left(\frac{1-\sqrt{1-2t}}{\mu}\right),\qquad t<1/2, ]

take logs:

[ \log M(t) = -\frac{1}{2}\log(1-2t) + \frac{1-\sqrt{1-2t}}{\mu}. ]

Differentiate:

[ \frac{d}{dt}\log M(t) = \frac{1}{1-2t} + \frac{1}{\mu\sqrt{1-2t}}. ]

Evaluating at (t=0) yields:

[ \mathbb{E}[X] = M’(0) = \left.\frac{d}{dt}\log M(t)\right|_{t=0} = 1 + \frac{1}{\mu}. ]

6.2 Variance#

Differentiate again:

[ \frac{d^2}{dt^2}\log M(t) = \frac{2}{(1-2t)^2} + \frac{1}{\mu}(1-2t)^{-3/2}. ]

At (t=0) this is the second cumulant (the variance):

[ \mathrm{Var}(X) = 2 + \frac{1}{\mu}. ]

6.3 Likelihood and MLE (standardized form)#

For i.i.d. data (x_1,\dots,x_n) with (x_i>0), the log-likelihood is

[ \ell(\mu) = \sum_{i=1}^n \log f(x_i\mid\mu) = -\frac{n}{2}\log(2\pi) - \frac{1}{2}\sum_i\log x_i - \sum_i\frac{(1-\mu x_i)^2}{2\mu^2 x_i}. ]

Expand the quadratic term:

[ \frac{(1-\mu x)^2}{2\mu^2 x} = \frac{1}{2\mu^2 x} - \frac{1}{\mu} + \frac{x}{2}. ]

So (dropping (\mu)-free constants) the objective becomes:

[ \ell(\mu) = \frac{n}{\mu} - \frac{1}{2\mu^2}\sum_{i=1}^n\frac{1}{x_i} + \text{const}. ]

Setting (\ell’(\mu)=0) yields a closed-form MLE:

[ \hat{\mu}{\mathrm{MLE}} = \frac{1}{n}\sum{i=1}^n \frac{1}{x_i}. ]

def recipinvgauss_loglikelihood(x: np.ndarray, mu: float) -> float:
    x = np.asarray(x, dtype=float)
    mu = float(mu)
    if mu <= 0 or np.any(x <= 0):
        return -np.inf
    return float(np.sum(recipinvgauss_logpdf(x, mu)))


def recipinvgauss_mle_mu(x: np.ndarray) -> float:
    x = np.asarray(x, dtype=float)
    if np.any(x <= 0):
        raise ValueError('all observations must be > 0')
    return float(np.mean(1.0 / x))


# Demonstrate MLE on synthetic data
mu_true = 1.7
x = stats.recipinvgauss(mu_true).rvs(size=5000, random_state=rng)
mu_hat = recipinvgauss_mle_mu(x)

mu_true, mu_hat

(1.7, 1.7103637479819758)

7) Sampling & Simulation (NumPy-only)#

A practical sampler uses the relationship:

[ X \sim \mathrm{RIG}(\mu) \iff \frac{1}{X} \sim \mathrm{InvGauss}(\mu). ]

So the plan is:

  1. sample (Y \sim \mathrm{InvGauss}(\mu)) using a NumPy-only algorithm

  2. return (X = 1/Y)

7.1 Inverse Gaussian sampler (Michael–Schucany–Haas)#

A widely used method for (Y\sim\mathrm{IG}(\mu,\lambda)) is:

  1. draw (V \sim \mathcal{N}(0,1)) and set (W = V^2)

  2. compute

[ Y^* = \mu + \frac{\mu^2 W}{2\lambda} - \frac{\mu}{2\lambda}\sqrt{4\mu\lambda W + \mu^2 W^2} ]

  1. draw (U\sim\mathrm{Uniform}(0,1)) and set

[ Y = \begin{cases} Y^*, & U \le \frac{\mu}{\mu+Y^*}\ \frac{\mu^2}{Y^*}, & \text{otherwise.} \end{cases} ]

SciPy’s standardized invgauss(mu) corresponds to (\lambda=1) in this notation.

def invgauss_rvs_numpy(mu: float, size: int, rng: np.random.Generator) -> np.ndarray:
    '''NumPy-only sampler for SciPy's standardized invgauss(mu).

    Uses the Michael–Schucany–Haas method with lambda=1.
    '''

    mu = float(mu)
    if mu <= 0:
        raise ValueError('mu must be > 0')
    if size <= 0:
        raise ValueError('size must be >= 1')

    v = rng.normal(size=size)
    w = v * v

    # lambda = 1
    y_star = mu + 0.5 * mu**2 * w - 0.5 * mu * np.sqrt(4.0 * mu * w + mu**2 * w**2)

    u = rng.random(size)
    y = np.where(u <= mu / (mu + y_star), y_star, mu**2 / y_star)
    return y


def recipinvgauss_rvs_numpy(mu: float, size: int, rng: np.random.Generator) -> np.ndarray:
    '''NumPy-only sampler for recipinvgauss(mu) via inversion of invgauss samples.'''

    y = invgauss_rvs_numpy(mu=mu, size=size, rng=rng)
    return 1.0 / y

# Validate sampler vs SciPy (quick KS test + moment check)
mu0 = 1.3
n = 50_000

samples_numpy = recipinvgauss_rvs_numpy(mu0, n, rng)
samples_scipy = stats.recipinvgauss(mu0).rvs(size=n, random_state=rng)

ks = stats.ks_2samp(samples_numpy, samples_scipy)

mc_mean = samples_numpy.mean()
mc_var = samples_numpy.var(ddof=0)

ks, (mc_mean, recipinvgauss_moments(mu0)['mean'], mc_var, recipinvgauss_moments(mu0)['var'])

(KstestResult(statistic=0.0027800000000000047, pvalue=0.9901116464472421, statistic_location=2.9112842485865573, statistic_sign=1),
 (1.7642504527853546,
  1.7692307692307692,
  2.761629618364975,
  2.769230769230769))

8) Visualization (PDF, CDF, Monte Carlo)#

We’ll visualize:

  • the PDF for multiple (\mu)

  • the CDF for multiple (\mu)

  • Monte Carlo samples vs the theoretical PDF and CDF

# PDF and CDF curves
mus = [0.5, 1.0, 2.5]
x = np.linspace(1e-4, 12, 700)

fig_pdf = go.Figure()
fig_cdf = go.Figure()

for mu in mus:
    fig_pdf.add_trace(go.Scatter(x=x, y=recipinvgauss_pdf(x, mu), mode='lines', name=f'μ={mu:g}'))
    fig_cdf.add_trace(go.Scatter(x=x, y=recipinvgauss_cdf(x, mu), mode='lines', name=f'μ={mu:g}'))

fig_pdf.update_layout(title='recipinvgauss PDF', xaxis_title='x', yaxis_title='f(x)')
fig_cdf.update_layout(title='recipinvgauss CDF', xaxis_title='x', yaxis_title='F(x)')

fig_pdf.show()
fig_cdf.show()

# Monte Carlo histogram vs PDF
mu0 = 1.3
n = 80_000
samples = recipinvgauss_rvs_numpy(mu0, n, rng)

x_grid = np.linspace(1e-4, np.quantile(samples, 0.995), 500)

fig = px.histogram(
    samples,
    nbins=70,
    histnorm='probability density',
    title=f'Monte Carlo samples vs PDF (n={n}, μ={mu0:g})',
    labels={'value': 'x'},
)
fig.add_trace(go.Scatter(x=x_grid, y=recipinvgauss_pdf(x_grid, mu0), mode='lines', name='true pdf'))
fig.update_layout(yaxis_title='density')
fig.show()

# Empirical CDF vs theoretical CDF
mu0 = 1.3
n = 30_000
samples = recipinvgauss_rvs_numpy(mu0, n, rng)

xs = np.sort(samples)
ys = np.arange(1, n + 1) / n

x_grid = np.linspace(1e-4, np.quantile(xs, 0.995), 500)

fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='empirical CDF'))
fig.add_trace(go.Scatter(x=x_grid, y=recipinvgauss_cdf(x_grid, mu0), mode='lines', name='true CDF'))
fig.update_layout(
    title=f'Empirical CDF vs true CDF (n={n}, μ={mu0:g})',
    xaxis_title='x',
    yaxis_title='CDF',
)
fig.show()

9) SciPy Integration (scipy.stats.recipinvgauss)#

SciPy exposes the distribution as scipy.stats.recipinvgauss with signature:

  • recipinvgauss(mu, loc=0, scale=1)

Recall the location–scale transform:

[ \texttt{pdf}(x;\mu,\mathrm{loc},\mathrm{scale}) = \frac{1}{\mathrm{scale}},f!\left(\frac{x-\mathrm{loc}}{\mathrm{scale}}\middle|\mu\right). ]

So:

  • changing scale multiplies the mean by scale and the variance by scale**2

  • changing loc shifts the support to (x > \mathrm{loc})

mu0 = 1.7
rv = stats.recipinvgauss(mu0)

x = np.linspace(1e-4, 6, 400)

pdf_vals = rv.pdf(x)
cdf_vals = rv.cdf(x)
samples = rv.rvs(size=5, random_state=rng)

pdf_vals[:3], cdf_vals[:3], samples

(array([0.    , 0.0001, 0.0132]),
 array([0.    , 0.    , 0.0001]),
 array([1.0811, 0.6593, 0.3056, 3.3773, 0.2918]))

# Fitting: compare SciPy's MLE (with fixed loc/scale) to the closed-form MLE
mu_true = 2.2
x = stats.recipinvgauss(mu_true).rvs(size=4000, random_state=rng)

mu_hat_closed = recipinvgauss_mle_mu(x)
mu_hat_fit, loc_hat, scale_hat = stats.recipinvgauss.fit(x, floc=0, fscale=1)

(mu_true, mu_hat_closed, mu_hat_fit, loc_hat, scale_hat)

(2.2, 2.2086896614467806, 2.208691406250003, 0, 1)

10) Statistical Use Cases#

10.1 Hypothesis testing (example)#

If you want to test a specific value (\mu=\mu_0), one option is a likelihood ratio (LR) test:

[ \Lambda = 2\left(\ell(\hat\mu) - \ell(\mu_0)\right). ]

Under regularity conditions, (\Lambda) is approximately (\chi^2_1) under the null.

10.2 Bayesian modeling (pattern)#

Because the log-likelihood simplifies to

[ \ell(\mu) = \frac{n}{\mu} - \frac{1}{2\mu^2}\sum_i\frac{1}{x_i} + \text{const}, ]

it is often convenient to reparameterize with (\theta = 1/\mu):

[ \ell(\theta) = n\theta - \frac{S}{2}\theta^2 + \text{const},\qquad S=\sum_i\frac{1}{x_i}. ]

As a function of (\theta), this is (up to constants) the log-density of a Gaussian (restricted to (\theta>0)). So a (truncated) normal prior on (\theta) yields a tractable posterior.

10.3 Generative modeling (idea)#

recipinvgauss is a positive distribution that can serve as a latent scale/precision. For example, sampling a latent (\tau>0) and then generating

[ Y\mid\tau \sim \mathcal{N}(0,\tau) ]

creates heavier-tailed marginals than a fixed-variance Gaussian.

# 10.1 Likelihood ratio test example
mu0 = 1.5
n = 800
x = stats.recipinvgauss(mu0).rvs(size=n, random_state=rng)

mu_hat = recipinvgauss_mle_mu(x)
ll_hat = recipinvgauss_loglikelihood(x, mu_hat)
ll_null = recipinvgauss_loglikelihood(x, mu0)

lr_stat = 2.0 * (ll_hat - ll_null)
# asymptotic p-value
p_value = stats.chi2(df=1).sf(lr_stat)

(mu0, mu_hat, lr_stat, p_value)

(1.5, 1.492655677690741, 0.012848466228206235, 0.9097522260451599)

# 10.2 Bayesian example on theta = 1/mu with a (truncated) Normal prior
# Prior: theta ~ Normal(m0, s0^2), restricted to theta>0.
# Likelihood in theta is proportional to exp(n*theta - (S/2)*theta^2).

x = stats.recipinvgauss(1.7).rvs(size=400, random_state=rng)

n = x.size
S = np.sum(1.0 / x)

# Likelihood corresponds to theta ~ Normal(mean=n/S, var=1/S) up to constants.
like_mean = n / S
like_var = 1.0 / S

m0, s0 = 0.6, 0.5  # prior mean/std on theta
prior_prec = 1.0 / (s0**2)
like_prec = 1.0 / like_var

post_prec = prior_prec + like_prec
post_mean = (prior_prec * m0 + like_prec * like_mean) / post_prec
post_std = np.sqrt(1.0 / post_prec)

# Sample theta from the (untruncated) Normal posterior and keep positive draws
n_draws = 50_000
theta_draws = rng.normal(loc=post_mean, scale=post_std, size=n_draws)
theta_draws = theta_draws[theta_draws > 0]
mu_draws = 1.0 / theta_draws

(mu_draws.mean(), np.quantile(mu_draws, [0.05, 0.5, 0.95]))

(1.619436549633346, array([1.4614, 1.6127, 1.7995]))

# 10.3 Generative modeling: a simple scale-mixture example
# Compare Y ~ Normal(0, 1) vs Y = sqrt(tau) * Z with tau ~ recipinvgauss(mu)

mu_tau = 2.0
n = 200_000

z = rng.normal(size=n)

# fixed-variance baseline
y_gauss = z

# scale-mixture
tau = stats.recipinvgauss(mu_tau).rvs(size=n, random_state=rng)
y_mix = np.sqrt(tau) * z

# compare tail quantiles
qs = [0.9, 0.95, 0.99, 0.995]
summary = {
    'q': qs,
    '|Y| Gaussian': [np.quantile(np.abs(y_gauss), q) for q in qs],
    '|Y| mixture': [np.quantile(np.abs(y_mix), q) for q in qs],
}
summary

{'q': [0.9, 0.95, 0.99, 0.995],
 '|Y| Gaussian': [1.6444642305320814,
  1.9596146024645902,
  2.5669629470726485,
  2.7937453680464364],
 '|Y| mixture': [1.9677723903719386,
  2.586557781439685,
  4.028211054677508,
  4.691210436430917]}

fig = go.Figure()

fig.add_trace(
    go.Histogram(
        x=y_gauss,
        nbinsx=200,
        histnorm='probability density',
        name='Normal(0,1)',
        opacity=0.6,
    )
)
fig.add_trace(
    go.Histogram(
        x=y_mix,
        nbinsx=200,
        histnorm='probability density',
        name=f'sqrt(tau)*Z, tau~RIG(mu={mu_tau:g})',
        opacity=0.6,
    )
)

fig.update_layout(
    title='Scale mixture increases tail mass',
    barmode='overlay',
    xaxis_title='y',
    yaxis_title='density',
)
fig.update_xaxes(range=[-6, 6])
fig.show()

11) Pitfalls#

  • Invalid parameters: mu must be (>0); scale must be (>0). Data must lie in the support.

  • Near-zero behavior: the PDF contains (\log x) and a (1/x) term in the exponent; evaluate via logpdf when possible.

  • CDF stability: the closed form uses differences of terms involving (\exp(2/\mu)); for extreme parameters use SciPy’s cdf, sf, logcdf, logsf for better numerical behavior.

  • loc/scale ambiguity: if you fit all parameters freely, location/scale can absorb structure; fix loc/scale when you know the data are standardized.

12) Summary#

  • recipinvgauss is a continuous distribution on (x>0) whose reciprocal is inverse Gaussian.

  • PDF: (f(x\mid\mu)=\frac{1}{\sqrt{2\pi x}}\exp\left(-\frac{(1-\mu x)^2}{2\mu^2 x}\right)).

  • Mean/variance: (\mathbb{E}[X]=1+1/\mu), (\mathrm{Var}(X)=2+1/\mu).

  • MGF exists for (t<1/2): ((1-2t)^{-1/2}\exp\left(\frac{1-\sqrt{1-2t}}{\mu}\right)).

  • Standardized MLE: (\hat\mu = \frac{1}{n}\sum_i 1/x_i).

  • Sampling is easy via NumPy-only inverse Gaussian sampling + reciprocal transform.

  • In practice, use SciPy’s scipy.stats.recipinvgauss for robust cdf/sf/fit and numerical entropy.

Contents